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A study of relative velocity statistics in Lagrangian perturbation 
theory with PINOCCHIO 
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ABSTRACT 

Subject of this paper is a detailed analysis of the PINOCCHIO algorithm for studying the 
relative velocity statistics of merging haloes in Lagrangian perturbation theory. Given a cos- 
mological background model, a power spectrum of fluctuations as well as a Gaussian linear 
density contrast field 5\ is generated on a cubic grid, which is then smoothed repeatedly with 
Gaussian filters. For each Lagrangian particle at position q and each smoothing radius R, the 
collapse time, the velocities and ellipsoidal truncation are computed using Lagrangian Pertur- 
bation Theory. The collapsed medium is then fragmented into isolated objects by an algorithm 
designed to mimic the accretion and merger events of hierarchical collapse. Directly after the 
fragmentation process the mass function, merger histories of haloes and the statistics of the 
relative velocities at merging are evaluated. We reimplemented the algorithm in C++, reco- 
vered the mass function and optimised the construction of halo merging histories. Comparing 
our results with the output of the Millennium simulation suggests that PINOCCHIO is well 
suited for studying relative velocities of merging haloes and is able to reproduce the pairwise 
velocity distribution. 

Key words: cosmology: large-scale structure, methods: analytical, numerical 



1 INTRODUCTION 

To account for structure formation one needs to develop techni- 
ques for studying the nonlinear evolution of perturbations. In the 
strongly nonlinear regime where the perturbation amplitudes ex- 
ceed unity (\6\ » 1), the linear approximation breaks down and 
has to be replaced by other approaches. To accentuate the need for 
this, one has to look at the interesting structures in the Universe, 
like galaxies or clusters of galaxies because they are highly non- 
linear. Since every theory of structure formation must be capable 
of describing the formation and evolution of non-linear objects, the 
major developments occurred in perturbation theory and numerical 
simulations. It has been understood that this aim simplifies when 
formulated in terms of Lagrangian coordinates rather than the stan- 
dard Eulerian ones, since the latter one relies on physical densities 
being small (Bouchet et al. 1992; Buchert 1992). So one assumes 
that the dynamics of gravitational clustering is described more sui- 
tably in terms of the displacement field D, which is in the Lagran- 
gian approach the only underlying fundamental field. The decisive 
difference to the Eulerian approach is that it is not based on the 
smallness of the density of the inhomogeneities and that one sear- 
ches for solutions of perturbed trajectories about the linear initial 
displacement Z) (1) (Bernardeau et al. 2002; Ehlers & Buchert 1997; 



e-mail: lavinia.heisenberg@unige.ch 



Sahni & Coles 1995). The fundamental point is that a small per- 
turbation of the Lagrangian particle paths carries a large amount of 
non-linear information about the corresponding Eulerian evolved 
observables, since the Lagrangian picture is intrinsically non-linear 
in the density field. Lagrangian perturbation theory reaches its limit 
of applicability when trajectories of particles cross and the mapping 
from the initial conditions to the evolved density field ceases to be 
unique. This defines the so-called orbit-crossing which is numeri- 
cally achieved by means of the ellipsoidal collapse approximati- 
on to the full Lagrangian perturbative expansion. These collapsed 
objects can then be grouped into disjoint haloes and their relative 
velocities can be measured. 

There is a number of applications of large simulated volumes 
with velocity information: For instance, simulations of redshift- 
space distortions in the galaxy correlation function, investigations 
of large-scale bulk flows (Zaroubi et al. 2002) which offer a possi- 
bility of testing modified gravity models (Ayaita et al. 2009), and 
testing the occurence of high-velocity merging events, such as the 
bullet-cluster, which is the topic of recent controversies (Hayashi & 
White 2006). Many cluster observables such as the strong-lensing 
cross-section (Fedeli et al. 2006; Torri et al. 2004) and the radio and 
X-ray luminosity (Randall et al. 2002) are boosted during the mer- 
ging process, and for treating these systems, statistical information 
about the merging configuration is needed. 

The cosmological model used is a spatially flat ACDM cos- 



2 L. Heisenberg, B.M. Schdfer and M. Bartelmann 



mology with Gaussian adiabatic initial perturbations in the cold 
dark matter (CDM) density field with the following cosmologi- 
cal parameter set (Q„o, S1 A , h, o" 8 ) = (0.25,0.75,0.73,0.9), iden- 
tical to that used in the Millennium simulation, to which we will 
compare our results. Our computational domain is a cubic box 
of size 256 Mpc/h with periodic boundary conditions filled with 
N 3 = 128 3 particles. 



2 COSMOLOGY 

In this first introductory section we summarize all the known for- 
mulas we use to compute the power spectrum and the Gaussian 
linear density contrast field. The linear CDM power spectrum P(k) 
describes the fluctuation amplitude of the Gaussian initial density 
field (5, (S(k)S(k')) = (2nfS D (k+k')P(k), and is given by the ansatz 



P(k) oc k"'T 2 (k), 



(1) 



with the transfer function T(k) which is approximated by (Bardeen 
et al. 1986), 

= Hl+aq) , + 2 + 3 + ( )4 r i (2) 
aq v ' 

In the transfer function, the wave number k = qT enters rescaled by 
the shape parameter T (Sugiyama 1995), 



T = Q. m h exp \-£l b 1 + 



V2A 



(3) 



For further treatment one introduces the smoothed density 
field Sr, which is averaged on the scale R, 



S R (k) = S(k)W R (k), 

and for the smoothed power spectrum 

P R (k) = \W R (k)\ 2 P(k). 



(4) 



(5) 



with the window function W R (r), satisfying f d 3 x W R (\x\) = 1. The- 
re exist several choices for the window function W R , a very com- 
mon one is the top-hat filter, which is given by 



W R (k) = 3 



sin (kR) - kR cos (kR) 
(kW) ' 



(6) 



Using the definitions introduced so far, it is now straightforward 
to calculate the variance <r 2 (R) = (6 2 r (x)) of the smoothed density 
field as 

, C d 3 k r d 3 k - , 

^ = TTj Pr( ® = \W R (k)\ 2 P(k). (7) 

J (2/r) 3 J (2n) 3 

which we use to fix the normalization to <t 2 (R = 8 Mpc//z) = <r 2 . 

3 LAGRANGIAN PERTURBATION THEORY 

In this section we summarize the idea and results of the Lagrangian 
perturbation theory, which we use to study the dynamical equati- 
ons of the considered non-linear system. The formation of CDM 
haloes involves highly non-linear dynamical processes which can 
only be followed in numerical simulations or by applying pertur- 
bative methods. As mentioned in the introduction, one important 
approach is Lagrangian Perturbation Theory. In order to describe 
the non-linear dynamical evolution of the particle trajectories up to 
third order Lagrangian coordinates the following three steps have 
to be followed: 



(i) Description of the mapping from the Eulerian to the Lagran- 
gian coordinates 

(ii) Transformation of the Eulerian fields to Lagrangian coordi- 
nates 

(iii) Application of perturbation theory to the Lagrangian equa- 
tions expressed by the displacement up to desired order 

The starting point is a pressureless, non-vortical, self-gravitating 
fluid with Newtonian gravity embedded in an expanding 
Friedmann-Lemaitre Robertson- Walker universe (Bouchet et al. 
1992). In the Lagrangian framework of fluid dynamics the relation 
between the Eulerian position x of a mass particle and the initial La- 
grangian position q is given by the displacement field (Zel'dovich 
1970): 



x(q, t) = q + D(q, t), 



(8) 



In the Lagrangian space the trajectories of the mass elements are 
fully described by the dynamical mappings x(q, t), starting from 
the initial positions q. 

There is a one-to-one correspondence between the Langrangi- 
an coordinate q and the Eulerian coordinate x for a cold non- 
collisional fluid, at least until the stage of caustic formation. Ex- 
pressed mathematically, this corresponds to the statement that the 
functional determinant J of the Jacobian of the mapping relation 
q — > x(q, t) is non-singular, 

7(3,0 = det(gUo, (9) 

which means that the mapping x(q, t) is invertible to q(x, t). It is 
evident that many particles coming from very different original 
positions will tend to arrive at the same Eulerian position during 
the highly non-linear evolution. As a consequence of that infinite- 
density regions (caustics) will form in Eulerian space. Hence the 
mapping from Lagrangian to Eulerian space becomes singular and 
the density infinite asp oc J~ l . Since the displacement D fully cha- 
racterizes the map between the Eulerian and the Lagrangian coor- 
dinates, the motion of the fluid elements are completely described 
in terms of it. One can now express the peculiar velocity, accele- 
ration and density contrast by inserting the displacement field D 
into the Euler and the continuity equations. Similarly the Eulerian 
irrotationality condition and the Poisson equation can be written in 
terms of the displacement (Catelan 1995) resulting in the dynamical 
equations for the displacements D. The fundamental question now 
is how to solve the dynamical equations for the displacements D. 
The irrotationality condition and the Lagrangian Poisson equation 
are exact equations in the Lagrangian description. It is undoubtedly 
very difficult to solve them in a rigorous way. A possible alternative 
is to seek for approximate solutions. The standard technique is to 
expand the trajectory D in a perturbative series, the leading term 
being the linear displacement which corresponds to the Zel'dovich 
approximation (Bouchet et al. 1992; Zel'dovich 1970). Approxima- 
ting the Lagrangian Poisson equation implies that the gravitational 
interaction among the particles of the fluid is described by the first 
few terms of a Taylor expansion. 

D(q, T) = gl (j) D m (q) + g 2 ( T ) D (2 \q) + ■■■ (10) 

D i '"\q) being the nth-order approximation. The dynamics of the 
evolution constrains in general both the temporal dependence as de- 
scribed by the functions g n , and the spatial displacements D {n \q). 
The displacement fields are computed considering the potential 
of a homogeneous ellipsoid in its principal axis frame ij/(q) = 
\{^iq\ + /l2<?2 + ^3^3) solving the Poisson and irrotationality condi- 
tion equtions in Lagrangian space, where the A t are the eigenvalues 



relative velocities with PINOCCHIO 3 



of the first-order deformation tensor >js M b(q)- The solutions can be 
found in (Catelan 1995) and (Monaco 1997a,b) to which we re- 
fer for all details. We use the results from (Monaco 1997a) for the 
calculation of collapse times J(q, g c ) = 



Tk + ¥4 AM6 >- 



*d)g 3 c =0. (11) 



Just in order to give some examples, the first order collapse time 
is given by g c ' 

(2) _ 7,<3+y7i 3 U 3 +6<i ) ) 

Monaco 1997a). 



4- whereas the second order solution results to 

■<3 



(and analogously for the third order, see 
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particle 




3.1 PINOCCHIO 

A compromise between simulations and analytical techniques is a 
perturbative approach describing the growth of haloes in a given 
numerical realization of a linear density field, such as the truncated 
Zel'dovich approximation and the PINOCCHIO algorithm (Taf- 
foni et al. 2002). PINOCCHIO (acronym for PINpointing Orbit- 
Crossing Collapsed Hierarchical Objects) is an algorithm for stu- 
dying the formation and evolution of dark matter haloes in a given 
initial linear density field, which we outline briefly in this section. 
It was first developed by Monaco et al. (2002). Local parameteri- 
sations to the dynamics are used to give precise predictions of the 
hierarchical formation of dark matter haloes when the correlations 
in the initial density field are properly taken into account. 

This modus operandi enables the automatic generation of a 
large ensemble of accurate halo merging histories and additionally 
delivers their spatial distribution. Likewise, the approach can be ef- 
ficiently applied for generating the input for galaxy formation mo- 
dels since the properties of the halo population are of fundamental 
importance for understanding galaxy formation and evolution. 

PINOCCHIO consists of two steps which determine the hier- 
archical formation of haloes through accretion and merging: 

• The first step handles the definition of the collapse time. Her- 
eby, orbit-crossing will be identified as the instant when a mass 
element undergoes collapse, without the need to introduce a free 
parameter. Orbit-crossing is numerically calculated by means of 
the ellipsoidal collapse approximation to the full Lagrangian per- 
turbative expansion, as discussed in the previous section (Monaco 
1995). 

• The second step groups the collapsed particles into disjoint 
haloes, applying an algorithm similar to that used to identify haloes 
in n-body simulations. 

As explained in the previous subsection, the density diverges as 
the Jacobian determinant vanishes, corresponding to the formation 
of a caustic, which means that particle trajectories cross and the 
transformation x — > q becomes multi-valued. Since the density be- 
comes very high at orbit crossing, this event will be identified as 
the collapse time. In this manner, collapse is very easy to compute 
using Lagrangian Perturbation Theory which remains valid up to 
that particular instant but breaks down afterwards. 

The hierarchical formation of objects is done due to the grou- 
ping of orbit-crossing particles into haloes by tracing the merging 
processes for each particle individually. Briefly, two main proces- 
ses contribute to the hierarchical clustering: The first one is the ac- 
cretion of particles onto haloes and the second one the merging of 
haloes. For this purpose, the particles of the realization are sorted 
in chronological order of collapse. Starting with the earliest collap- 
se time the particles are assigned either to a halo or to filaments at 
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Figure 1. Cases of the fragmentation process: The top panel shows the six 
Lagrangian neighbours of a given particle, the second panel illustrates how 
this particle accretes onto a neighbouring halo, the third panel depicts the 
merging of two haloes and successive accretion and, if there is no accretion, 
the particle is marked as belonging to a filament in the bottom panel. 
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the corresponding collapse times. In order for a collapsing partic- 
le to be accreted by a halo, it must fulfill a number of conditions. 
One of them is that the candidate halo must already contain one of 
its 6 nearest neighbours in the Lagrangian space of initial condi- 
tions. The fragmentation process contains 4 different cases which 
are sufficient for performing the identification and merging history 
of haloes (Fig. 1). 

If none of the 6 Lagrangian neighbours have collapsed, then 
the particle is a local maximum of the inverse collapse time. This 
particle is a seed for a new halo having the unit mass of the particle 
and is created at the particle's position. Obviously, the earliest par- 
ticle to collapse is the first halo. 

In the case that the collapsing particle touches only one halo, then 
the accretion condition, if the halo is close enough, is checked. 
When the accretion condition is satisfied, then the particle is ad- 
ded to the halo, otherwise it is marked as belonging to a filament. 
Thus, not every collapsed particle becomes a relaxed halo but can 
be tagged as a filament. The particles that only touch filaments are 
marked as filaments as well. If the collapsing particle has more than 
one touching halo as Lagrangian neighbour, then the merging con- 
dition is checked for all halo pairs. Pairs that satisfy the conditions 
are merged. The accretion condition for the particle is checked for 
all touching haloes both before and after merging. In the case that 
the particle can accrete to both haloes, without the haloes merging, 
it accretes onto that halo for which the distance d in units of halo 
size R N is smaller. It may happen that particles fail to accrete even 
though the haloes merge. 

If the collapsing particle does not accrete onto the candidate ha- 
loes in the case they are too far, it becomes a filament. But later 
for this filament particle there is still the possibility to accrete when 
its neighbour particle accretes onto a halo. This is done in order to 
mimic the accretion of filaments onto the haloes. Notice that up to 
5 filament particles can flow into a halo at each accretion event. 

3.2 Velocity statistics 

In the Eulerian specification of the flow field, the flow quantities 
are depicted as functions of fixed position q and time t. The pecu- 
liar velocity expressed in terms of the displacement field D in the 
Lagrangian specification of the flow field would then be 



V(t) = i(f) = g(t)D(q) 



(12) 



Now let us study the time evolution of two elements positioned at 
the Lagrangian coordinates q a and q b . The Eulerian positions and 
peculiar velocities at these two points are then given as 

x a = X a (q a ,t) = q a + g(t)D(q a ) = q a + g(t)D a , (13) 

x b = x b (q b , t) = q b + g(t)D(q b ) = q b + g(t)D b , (14) 

v a (f) = v a (q a ,t) = g{t)D a , (15) 

v b (t) = v b (q b ,t) = g(t)D b , (16) 

The relative velocity of the two particles is given by 

V ab = V b (t) - V a (f) = g(t)(D b - D a ) = W||(f) + v ± (t), (17) 

where t/||(f) and v ± (t) stands for the components parallel and per- 
pendicular to r ab (t) = x b (t) - x a (t). We can now compute the pro- 
bability distribution function (PDF) of the pairwise peculiar velo- 
city v with separation s from the initial PDF as (Seto & Yokoyama 
1998): 

P(v, s, t) = -^-j J 4nr 2 dr dv lv dv ±:si du ±yi p(u llh u ±li , v Xyi ; r) 

xS( S -r ab (t))6(v-v lv (t)) (18) 




Figure 3. the visualisation of the configuration 

where p(v^, v ±xi , v ±yi ; r) is the initial PDF which depends only on 
V\\i and v n = yjv 2 ±tj + u 2 ±yj where u ±yi and u ±yi are the two com- 
ponents of v ±j perpendicular to each other. The subscript i denotes 
quantities at some initial time. 

Putting the Eulerian positions and peculiar velocities into the 
above PDF gives: 

P(v, s , t) oc J" rdr p(v*p V*f, r), 

where u* and i/*. stand for 
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rg 

.2 



5' 1 
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g 
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(19) 

(20) 
(21) 
(22) 



We have to specify the initial PDF p(vt, v* } ; r) in order to be able 
to compute the desired PDF P(v,s,t). As we were dealing with 
the longitudinal mode, the peculiar velocity in the linear regime is 
related to the 6(x, t) and accordingly to its Fourier transform 6k(t) 



v(x, t) = i 



li f A 

giJ k> 



(2nf ■ 



(23) 



It is important to emphasize that the initial pairwise peculiar velo- 
cities are Gaussian distributed like the initial density fluctuations. 
From the velocity correlation tensor {v a v b } one obtains the projec- 
tions: 



< ;, ll !/ ll) 



E 

u.b 



(24) 



(25) 



Thus, the two-point correlation functions are given by (Gorski 
1988): 

\2 



(VxiVxi) 



± (t) J dkP,(k){l - 3j (kr) + 6^) (26) 



(t) 1-3 



Mkr) 
kr 



Finally, the initial probability distribution function is given by 



p(v*, v* ±i ; r) 



T = 



+ 



2Fn(r) 2YAr) 



(27) 



(28) 



where (v^v^) = Y^(r) and (u ±i v ±i ) = Y ± (r). Therefore, we can now 
obtain the desired PDF P(v, s, t) through the integration P(v, s, t) oc 
rdr p(v*^, r). We relate the distance r in the analytical mo- 
del to the distance of two haloes at the time of merging. Furthermo- 
re we substitute for the power spectrum a filtered spectrum smoo- 
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Figure 2. Flow chart of the code: The main 3 blocks of the code are shown. For a given set of parameters (centre) we compute the collapse times (top left) and 
for each collapsing particle we apply the fragmentation procedure (bottom). Finally we analyse the statistical properties of the fragmented objects (top right). 



thed at a wavelength which corresponds to the halo mass. Compu- 
ting the velocity distribution at a distance s corresponding to the 
merging condition yields distributions very similar to those obtai- 
ned by PINOCCHIO. In particular we confirm the trend of steeper 
distributions at lower mass ratio. 

As expected the velocity distribution for the analytical PDF is 
shallower than the distribution from PINOCCHIO data. This re- 
flects the fact that merging processes in PINOCCHIO conserve 
momentum but not energy since merging is an inelastic collisi- 
on. At each merging the velocity of the final halo is given by 
Vf = {vinii + V2in2)/(mi + mi) since the momentum is conser- 
ved. Because of the energy loss high velocities do not appear in the 
PINOCCHIO velocity distribution and therefore the curve decrea- 
ses faster. But the loss of energy does not occur in the analytical 
configuration and therefore the curve is shallower. 



4 RESULTS 

In the following we present the results from a simulation for 128 3 
particles in a box of side length 256 Mpc//7, carried out in the fra- 
mework of Lagrangian perturbation theory with our C++ reimple- 
mentation of the PINOCCHIO code for following the merging ac- 
tivity of the large-scale structure. There is always the problem of 
the coherence of the velocity field on large scales. For that rea- 
son, large simulation boxes are always preferrable for the investi- 
gation of velocity statistics. Due to nonlinear mode-coupling even 
large modes influence the dynamics on small scales (Cole 1997), 
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Figure 4. The analytical pair velocity probability density for three intervals 
in mass ratio: low mass ratio 0.0 < log(M)/M2) < 0.7 (red line), medium 
mass ratio 0.7 < log(Mi/M2) < 1.3 (green line) and high mass ratio 1.3 < 
log(Mi/M 2 ) < 2.0 (blue line). 



although this coupling decreases proportional to the inverse wave- 
number or faster in classical perturbation theory. This difficulty can 
in principle be treated with PINOCCHIO using the method descri- 
bed by (Monaco et al. 2005) which consists in adding large-scale 
modes to the Gaussian initial conditions, furthermore a volume of 
size 256 Mpc//r has been demonstrated to sample the velocity field 
sufficiently well, as velocity differences are much less dependent on 
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large-scale modes than velocities: Heuristically, the velocity diffe- 
rence statistic is closely related to that of the velocity field diver- 
gence, and has therefore the same properties as the density field 
itself, which is well sampled by our volume. 

The first step was to create a realisation of a Gaussian densi- 
ty field for a specified ACDM spectrum. From the density contrast 
the gravitational potential and Zel'dovich tensor were derived. Dif- 
ferentiations were performed again in Fourier space allowing us 
to recover the quantities with minimum noise. For each point q of 
the Lagrangian initial coordinates and each smoothing radius R, 
the collapse time, determined by the time at which the particle is 
predicted to enter a high-density multi-stream region, and ellipsoi- 
dal truncation were then computed using Lagrangian perturbation 
theory. Determining the eigenvalues of the Zel'dovich tensor, it is 
straightforward to calculate the collapse times evaluating the ex- 
pressions of the individual orders (11). For each particle only the 
earliest collapse time is recorded with the corresponding smoothing 
radius and velocity. The third order Lagrangian prescription impro- 
ves the collapse times. The first axis to collapse is the one corre- 
sponding to the smallest A 3 eigenvalue indicating the convergence 
and the fact that the Zel'dovich approximation makes the largest 
contribution to collapse dynamics. The collapsed medium is then 
fragmented into isolated objects using the algorithm which we des- 
cribed above to mimic the accretion and merger events of hierarchi- 
cal collapse. We distinguish between collapsed particles belonging 
to relaxed haloes or to lower-density filaments using the accretion 
and merging conditions. Once the fragmentation process is com- 
pleted we studied the relative velocities of the haloes at merging. 

Figs. 7 and 8 illustrate structure formation dynamics in the 
Lagrangian picture: Matter is transported out of initially underden- 
se regions and is accumulated in superstructures, where merging is 
prominent due to the high particle density. In fact, the most mas- 
sive objects are found in regions of converging velocities, a nice 
example of which can be found in the upper right front corner of 
the simulation cube. One sees coherent flow patterns on the scale 
of the correlation length of the density field. In the vicinity of mas- 
sive structures one can observe larger relative velocities compared 
to underdense regions. This can be traced back to the fact that the 
velocity field of an overdense region has a larger variance compa- 
red to the cosmological average. The statistics of the velocity field 
is translated to that of the haloes by imposing momentum conser- 
vation in the merging process. Therefore, the figure confirms the 
expectation of high pairwise velocities in overdense regions. 

In order to compare our results with the Millennium Simu- 
lation (Springel et al. 2005) we have applied the same merging 
conditions to the Millennium data. In Figure 5 we compare the 
probability distribution of the relative velocity at the time of mer- 
ging, for three different mass ratio intervals: approximately equal 
masses log(M t /M 2 ) < 0.7, intermediate values for the mass ratio 
0.7 < log(Mi/M 2 ) < 1.3 and high mass ratios log(Mi/M 2 ) > 1.3. 
The data is further split into three redshifts: z = 0, 1,2, while at 
higher redshifts the numbers of massive haloes is not sufficient for 
deriving the probability density as statistical error bars are too large 
to draw conclusions. The subdivision of the velocity range between 
0-1000 km/s/h allows us to investigate the features of the velocity 
distribution while obtaining reasonable statistical error bars, and 
we normalise all histograms to unity. The error bars of the PINOC- 
CHIO runs (V = (256 Mpc/h) 3 ) are small enough to allow quantita- 
tive conclusions and are comparable to those from the Millennium- 
data (V = (256 Mpc/h) 3 ). For drawing velocity distributions we 
selected halos from mass intervals of identical width ranging from 
4 x lQ n M e /h to 1 x 10 14 MJh. 



mass function 




Figure 6. comparison of the PINOCCHIO mass function with the analytical 
Press-Schechter, Sheth-Tormen and Jenkins mass functions 



For controlling the time-stepping, GADGET (which was used 
to carry out the Millennium simulations) has a feature which syn- 
chronizes the individual time-steps for each particle prior to writing 
a simulation output at a pre-specified redshift. The Pinocchio simu- 
lation was stopped at the identical redshift, and given this approach, 
we do not expect difficulties to arise due to the time-discreteness of 
the Millennium output. Furthermore, as long as the structures are 
in the limit, where perturbation theory is applicable, which means 
that the trajectories of the particles can be extrapolated from the 
initial conditions, the choice of the time-stepping should not mat- 
ter much, in contrast to nonlinear structures with large gradients 
perpendicular to the particle trajectory. 

A Metropolis-Hastings algorithm was used for checking the 
dependence of the mass function on the PINOCCHIO-parameters 
and we were able to reproduce mass functions in agreement with 
analytical functions for the proposed set of parameters. The influ- 
ence of the parameter choice on the shape of the mass function 
was verified in our implementation, and the total number of haloes 
formed corresponds well to analytical expectations and the Mill- 
ennium simulation. For most of the tests presented in this paper 
we used a Opteron 285 2.4 GHz 64 bit computer with a total me- 
mory of 8GB. For example, resampling the initial conditions onto 
a 256 3 grid, the first step of initialising the Gaussian random field 
and computing the orbit-crossing requires about 30 minutes of CPU 
time and the second part of identifying the haloes takes only a few 
additional minutes. Fig. 6 illustrates the agreement of the PINOC- 
CHIO mass function with the analytical expectation. 

Quite generally, the figures reproduce the basic behaviour ex- 
pected from analytical arguments as outlined in Sect. 3.2. The dis- 
tributions are steeper at lower mass ratio in both the Millennium 
data set as well as in PINOCCHIO. The shape of the distributions 
between PINOCCHIO and the Millennium simulation is very simi- 
lar, particularly well at low redshifts, with the curves for the inter- 
mediate and the high mass ratio peaking at almost identical values 
for the relative velocity. At high redshifts PINOCCHIO underesti- 
mates the velocities in medium and high mass ratio mergers by a 
small amount. This underestimate is a known feature of Lagrangian 
codes and is discussed in the appendix of Monaco et al. (2005). All 
curves terminate at velocities of ^ lOOOkm/s underlining the spar- 
city of high-velocity mergers (Hayashi & White 2006). Naturally 
we expect the distributions of the Millennium simulation to drop 
faster than corresponding distributions from PINOCCHIO, because 
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Figure 5. Probability distributions of velocities conditional on halo mass ratio for redshift z = (first row), z = 1 (second row) and z = 2 (third row), comparing 
PINOCCHIO (left column) with the Millennium simulation (right column). 



the latter treats merging processes as an inelastic collision and does 
not follow the dissipative dynamics inside haloes. Given the good 
agreement, we believe that halo velocity catalogues for a number 
of applications such as redshift space distortions, large-scale bulk 



flows and merging processes can be reliably derived from Lagran- 
gian codes. 
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Figure 8. Similar to Fig. 7 only without isodensity contours for making the velocity field more visible. Logarithmic halo mass is indicated by size, and the 
arrow length is proportional to the halo velocity. 



relative velocities with PINOCCHIO 9 



5 SUMMARY 

The aim of this paper is an investigation of merging processes in 
the cosmic large-scale structure in the framework of Lagrangian 
perturbation theory and to compare the results for the pairwise ve- 
locity distribution obtained with an adaptation of the PINOCCHIO 
algorithm to n-body data. 

• It is comparatively simple in PINOCCHIO to construct merger 
trees as it is not necessary to identify haloes in the particle distribu- 
tion and as one can directly follow the merging processes between 
haloes, such that the haloes in PINOCCHIO correspond to friends- 
of-friends particle groups in n-body simulations. 

• Halo properties, such as the distribution of masses, spins and 
now also merger trees and pairwise velocities can be reliably deri- 
ved from a Lagrangian code, with the differences in the distribution 
being smaller than the statistical error bars. Hereby, we have inves- 
tigated the derivation of the mass function from PINOCCHIO and 
used similiar parameters as quoted by Monaco et al. (2002) for our 
studies of the velocity distributions. We have also investigated that 
the velocities show the correct scaling with cosmological parame- 
ters, and the correct scaling with halo mass ratio in comparison to 
an analytical calculation. 

• In comparison to an analytical pairwise velocity distribution, 
PINOCCHIO is able to reproduce the trend of shallower distribu- 
tions with increasing halo mass ratio. At high velocities, however, 
PINOCCHIO exhibits a steeper behaviour compared to that predic- 
ted by the analytical calculation at fixed mass ratio, which is explai- 
ned by the fact that merging processes in PINOCCHIO are treated 
as inelastic collisions with conserved momentum but not conserved 
energy. Because of this energy loss, high velocities are not present 
in PINOCCHIO data and the distribution is steeper. 

• We find a general agreement between the velocity distributi- 
ons of PINOCCHIO and the Millennium simulation, both in terms 
of relative numbers and values for the absolute velocity. Additional- 
ly, the scaling with redshift and mass ratio between merging haloes 
behaves very similarly. The peaks of the distributions at low reds- 
hifts coincide with each other, and although distributions from the 
Millennium simulation terminate earlier, this feature is not unex- 
pected as PINOCCHIO does not treat the dissipative dynamics of 
merging haloes. Again, the correct dependence of pairwise velocity 
with halo mass ratio is recovered. 

• PINOCCHIO, relying on a phenomenological description of 
the merging process of two haloes and combining the individual 
momenta in an inelastic collision is very fast compared to rc-body 
codes, which allows sweeps in the parameter space relevant to pe- 
culiar velocities, i.e. Cl m , <x 8 and the dark energy parameters, for 
which we have verified the basic relations expected by linear struc- 
ture formation. 

Further questions include the environment-dependence of velocity 
statistics, which is comparatively easy to do in Lagrangian pertur- 
bation theory. A very useful discriminant for this purpose is the 
number of positive eigenvalues of the shear tensor. One would 
expect smaller velocities inside voids and larger velocities in su- 
percluster regions. In fact, this dependence can already be seen in 
Fig. 7. Other extensions include the investigation of two-point sta- 
tistics of the velocity field, and to answer questions related to ve- 
locity statistics (Regos & Szalay 1995). In a future paper, we will 
investigate the dependence of the strong lensing cross-section of 
merging systems with velocities drawn from PINOCCHIO simu- 
lated volumes, and its dependence on the choice of cosmological 
parameters. 



Anaglyphic versions of Figs. 7 and 8 are available on request 
from the authors. 
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